Allyson Terry’s study on immunologic response to liver allograft injury

Statement of statistical hypothesis(es):

  1. What cytokines are significantly different between IRI- and IRI+ patients in the larger cohort of 100 liver transplant patients?
  1. Do these cytokines influence IRI individually?
  2. Do these cytokines influence IRI together?
  3. Are any of these cytokines significantly different at more than one timepoint?

To determine if cytokines influence IRI individually, binary logistic regression will be conducted with IRI injury as the outcome variable, with difference between timepoints in cytokine concentrations as covariates.

To address what cytokines influence IRI together, a correlation plot was constructed to visualize related cytokines.

To address whether cytokines significantly differ at more than one timepoint, separate one-way anovas with visit type as the factor variable will be conducted with each cytokine as the outcome variable. Furthermore, two-way anovas will be conducted for each cytokine concentration with visit type and IRI status as factor variables to explore if time and IRI status explain variation in cytokine concentration.

  1. What cytokines in the pre-transplant blood sample specifically correlate with:
  1. Immune cell infiltration after transplant? (Inflammation score)
  2. Dead cells in the liver after transplant? (Necrosis score)

Background information:

Ischemia-reperfusion injury (IRI) is when an organ loses and then regains access to blood flow, commonly occurring with organ transplantation. The immune response elicited by IRI can be severe enough to negatively affect the function and survival of the donated organ. One of two main ways that the cells in the body communicate is through cytokines, which are secreted proteins that are present in the blood of transplant patients and can influence the degree of inflammation that occurs during and after IRI. Two characteristics that doctors use to determine if a patient is IRI+ or IRI- are the inflammation score, which is the amount of immune cells in the liver after transplant and the necorsis score, which is the amount of dead cells in the liver after transplant. These characteristics of inflammation and necrosis are scored from 0 to 4, with 0 being the least severe and 4 being the most severe.

Sample data description:

The study includes 100 liver transplant patients, 42 of which are IRI+ and 58 IRI-. The analysis includes 23 cytokines. The data set also includes the inflammation score, necrosis score, and visit type for each patient.

Variable assessment:

data <- readxl::read_xlsx("./20210521_402a.xlsx", col_names = T, skip = 1)

str(data)
## tibble [300 × 29] (S3: tbl_df/tbl/data.frame)
##  $ Sample #       : num [1:300] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Patient ID     : num [1:300] 4 4 4 5 5 5 6 6 6 7 ...
##  $ IRI Score (–,+): num [1:300] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Inflammation   : num [1:300] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Necrosis       : num [1:300] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Visit Type     : chr [1:300] "PO" "LF" "W1" "PO" ...
##  $ CYTO1          : chr [1:300] "3.62" "8.57" "70.98" "13.75" ...
##  $ CYTO2          : chr [1:300] "18.399999999999999" "361.15" "162.97" "120.96" ...
##  $ CYTO3          : num [1:300] 293 204 206 755 737 ...
##  $ CYTO4          : chr [1:300] "70.069999999999993" "68.459999999999994" "12.56" "2.63" ...
##  $ CYTO5          : chr [1:300] "<0.54↓" "2.76" "3.6" "7.51" ...
##  $ CYTO6          : num [1:300] 101 260 1558 364 1079 ...
##  $ CYTO7          : num [1:300] 22.1 3065 50.1 44.3 240.4 ...
##  $ CYTO8          : chr [1:300] "<2.26↓" "3.21" "<2.26↓" "8.68" ...
##  $ CYTO9          : chr [1:300] "7.41" "5.75" "20.81" "64.84" ...
##  $ CYTO10         : chr [1:300] "1.99" "10.62" "9.67" "8.9700000000000006" ...
##  $ CYTO11         : chr [1:300] "58.41" "187.22" "1478" "149.26" ...
##  $ CYTO12         : chr [1:300] "1.56" "3.06" "1.19" "8.74" ...
##  $ CYTO13         : chr [1:300] "25.57" "1210" "39.07" "127.7" ...
##  $ CYTO14         : chr [1:300] "41.76" "29.2" "97.12" "169.95" ...
##  $ CYTO15         : chr [1:300] "<1.45↓" "1.9" "1.98" "5.51" ...
##  $ CYTO16         : chr [1:300] "<1.09↓" "1.39" "<1.09↓" "<1.09↓" ...
##  $ CYTO17         : chr [1:300] "<1.38↓" "<1.38↓" "<1.38↓" "<1.38↓" ...
##  $ CYTO18         : chr [1:300] "261.95" "107.07" "671.93" "2177" ...
##  $ CYTO19         : chr [1:300] "5.0599999999999996" "3.49" "15.65" "32.85" ...
##  $ CYTO20         : chr [1:300] "71.97" "298.99" "98.69" "203.47" ...
##  $ CYTO21         : chr [1:300] "<2.07↓" "13.49" "6.1" "14.9" ...
##  $ CYTO22         : num [1:300] 37.8 38.1 43.2 103.3 61.1 ...
##  $ CYTO23         : num [1:300] 17.66 18.43 15.37 22.27 2.67 ...
dataclean <- data %>% 
  rename(sample = "Sample #",
         patientid = "Patient ID",
         IRI = "IRI Score (–,+)",
         visittype = "Visit Type") %>% 
  mutate(across(starts_with("CYTO"),
                ~ dplyr::case_when(
                  str_starts(.x, "<") == T ~ as.numeric(str_sub(.x, 2, -2)),
                  str_starts(.x, ">") == T ~ as.numeric(str_sub(.x, 2, -2)) + 1,
                  str_starts(.x, "1|2|3|4|5|6|7|8|9|0") == T ~ as.numeric(.x)
                )))
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.
## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following object is masked from 'package:purrr':
## 
##     transpose
isidata <- fread(file = "dataclean.csv", header = T, stringsAsFactors = F)
(isi_tble <- tibble(isidata))

Create patient characteristics table.

#namess
library(gtsummary)

CYTO <- names(isi_tble[7:29])  #vector of cytokine variable names 

demographic_table <- isi_tble %>% 
   mutate(IRI = ifelse( IRI == 0, "IRI-", "IRI+")) %>%
  select(-patientid, -sample) %>% 
  tbl_summary(by = IRI, 
              statistic = CYTO ~  
              "{median} ({p25}, {p75})",
              digits = all_continuous() ~ 2,
              label = list(visittype ~ "Visit Type"),
              missing_text = "Missing") %>% 
  add_overall() %>% 
  bold_labels() %>%
  add_p() %>% 
  bold_p() %>% 
  modify_header(label ~ "**Table 1. Patient Characteristics**" )
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(CYTO)` instead of `CYTO` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
demographic_table
Table 1. Patient Characteristics Overall, N = 3001 IRI-, N = 1741 IRI+, N = 1261 p-value2
Inflammation <0.001
0 27 (9.0%) 24 (14%) 3 (2.4%)
1 177 (59%) 126 (72%) 51 (40%)
2 78 (26%) 21 (12%) 57 (45%)
3 18 (6.0%) 3 (1.7%) 15 (12%)
Necrosis <0.001
0 75 (25%) 72 (41%) 3 (2.4%)
1 111 (37%) 99 (57%) 12 (9.5%)
2 72 (24%) 3 (1.7%) 69 (55%)
3 36 (12%) 0 (0%) 36 (29%)
4 6 (2.0%) 0 (0%) 6 (4.8%)
Visit Type >0.9
LF 100 (33%) 58 (33%) 42 (33%)
PO 100 (33%) 58 (33%) 42 (33%)
W1 100 (33%) 58 (33%) 42 (33%)
CYTO1 8.70 (4.30, 15.69) 8.94 (4.04, 15.62) 8.62 (4.83, 15.55) 0.6
CYTO2 128.04 (66.54, 359.61) 122.28 (66.32, 360.64) 129.10 (69.68, 352.29) 0.9
CYTO3 229.10 (151.17, 324.17) 232.78 (151.52, 326.69) 221.18 (149.24, 312.09) 0.7
CYTO4 24.62 (11.54, 45.87) 23.44 (13.08, 43.24) 25.07 (9.51, 47.88) 0.9
CYTO5 3.79 (1.53, 8.85) 3.79 (1.76, 8.32) 3.80 (1.33, 9.34) 0.8
CYTO6 271.99 (152.17, 459.83) 277.15 (172.53, 509.27) 262.06 (140.55, 397.65) 0.2
CYTO7 88.35 (44.57, 243.92) 94.28 (44.34, 276.30) 84.80 (46.80, 198.98) 0.7
CYTO8 3.33 (2.01, 7.73) 3.36 (2.00, 7.79) 3.29 (2.10, 7.70) 0.8
CYTO9 15.16 (5.07, 46.98) 15.32 (4.51, 48.28) 14.98 (6.52, 42.99) 0.8
CYTO10 10.12 (6.03, 15.21) 9.62 (5.75, 14.88) 10.27 (6.60, 16.62) 0.3
CYTO11 92.06 (44.36, 178.51) 101.12 (46.94, 189.43) 68.12 (37.17, 158.68) 0.037
CYTO12 2.86 (1.58, 6.04) 2.66 (1.58, 6.16) 2.99 (1.56, 5.90) 0.5
CYTO13 149.50 (69.78, 1,568.25) 145.42 (59.26, 1,468.50) 168.98 (82.33, 2,155.00) 0.3
CYTO14 66.06 (20.62, 189.05) 63.72 (17.14, 197.00) 71.24 (26.99, 179.49) 0.5
CYTO15 3.16 (1.51, 7.11) 3.00 (1.48, 7.26) 3.27 (1.67, 6.22) 0.8
CYTO16 1.10 (0.94, 1.60) 1.10 (0.81, 1.55) 1.10 (0.94, 1.92) 0.064
CYTO17 1.38 (1.33, 1.41) 1.39 (1.33, 1.41) 1.34 (1.06, 1.41) 0.018
CYTO18 431.17 (131.20, 1,351.75) 388.58 (109.13, 1,363.25) 491.25 (189.98, 1,316.75) 0.3
CYTO19 8.79 (3.79, 24.11) 7.98 (3.34, 25.27) 10.02 (4.52, 21.50) 0.5
CYTO20 152.46 (68.05, 306.70) 137.80 (61.75, 298.48) 163.06 (70.80, 336.67) 0.3
CYTO21 5.24 (2.79, 12.36) 5.24 (2.59, 13.66) 5.26 (3.28, 11.09) 0.8
CYTO22 58.27 (35.16, 107.86) 49.77 (31.09, 92.57) 75.42 (37.61, 117.85) 0.005
CYTO23 20.20 (13.20, 30.07) 18.40 (12.97, 28.57) 21.86 (13.45, 32.78) 0.065

1 n (%); Median (IQR)

2 Pearson's Chi-squared test; Fisher's exact test; Wilcoxon rank sum test

It should be noted that a significant difference between IRI- and IRI+ patients was detected for necrosis score (p <.001).

The significance test for cytokine concentration was conducted using a nonparametric method, so it is more robust to assumptions of normality. CYTO11, CYTO17, and CYTO22 were found to be significantly different. Next we visualize the distribution of cytokine concentration.

isi_tble %>% 
  select(CYTO1:CYTO6) %>% 
  gather() %>% 
  ggplot(aes(value)) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

isi_tble %>% 
  select(CYTO7:CYTO12) %>% 
  gather() %>% 
  ggplot(aes(value)) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

isi_tble %>% 
  select(CYTO13:CYTO18) %>% 
  gather() %>% 
  ggplot(aes(value)) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

isi_tble %>% 
  select(CYTO19:CYTO23) %>% 
  gather() %>% 
  ggplot(aes(value)) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## my attempt to facet in one plot. Does not look good
isi_tble %>% select(starts_with("CYT")) %>%
  gather() %>%
  ggplot(aes(value)) +
facet_wrap(~ key, scales ="free") +
# facet_grid_paginate(~ key, ncol = 3, nrow = 3, page = 2)
  geom_boxplot()

The Cytokine distributions are characterized by heavy right skew.

Now log transform all Cytokines

isi_tble %>% 
  select(CYTO1:CYTO6) %>% 
  gather() %>% 
  ggplot(aes(log(value))) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_histogram() + 
  labs(title = "Log Transformed") +
  xlab("log pg/mL")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

isi_tble %>% 
  select(CYTO7:CYTO12) %>% 
  gather() %>% 
  ggplot(aes(log(value))) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_histogram() + 
  labs(title = "Log Transformed") +
  xlab("log pg/mL")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

isi_tble %>% 
  select(CYTO13:CYTO18) %>% 
  gather() %>% 
  ggplot(aes(log(value))) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_histogram() + 
  labs(title = "Log Transformed") +
  xlab("log pg/mL")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

isi_tble %>% 
  select(CYTO19:CYTO23) %>% 
  gather() %>% 
  ggplot(aes(log(value))) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_histogram() + 
  labs(title = "Log Transformed") +
  xlab("log pg/mL")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

isi_tble %>% 
  select(CYTO1:CYTO6) %>% 
  gather() %>% 
  ggplot(aes(log(value))) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_boxplot() + 
  labs(title = "Log Transformed") +
  xlab("log pg/mL")

isi_tble %>% 
  select(CYTO7:CYTO12) %>% 
  gather() %>% 
  ggplot(aes(log(value))) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_boxplot() + 
  labs(title = "Log Transformed") +
  xlab("log pg/mL")

isi_tble %>% 
  select(CYTO13:CYTO18) %>% 
  gather() %>% 
  ggplot(aes(log(value))) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_boxplot() + 
  labs(title = "Log Transformed") +
  xlab("log pg/mL")

isi_tble %>% 
  select(CYTO19:CYTO23) %>% 
  gather() %>% 
  ggplot(aes(log(value))) + 
 facet_wrap(~ key, nrow = 3, scales = "free") + 
  geom_boxplot() + 
  labs(title = "Log Transformed") +
  xlab("log pg/mL")

  • Among the log transformed cytokine values, it is seen that that overall normality is improved.
  • CYTO12 and CYTO 19 appears to have a similar distribution when looking at the histogram, with large left skew and lack of density left of the peak density.
  • CYTO13 appears as if it may have another peak in the right tail. It is right skewed.
  • CYTO 15 appears to have several outliers in the higher range.
  • CYTO16 appears right skewed.
  • CYTO18 appears to have lower outliers.
  • CYTO17 has many values below level of detection, therefore the transformation did not affect overall normality.
  • CYTO21 appears to have a large right skew, though it should be noted that many values were below level of detection in CYTO 21.

It should be noted that due to many observations below detection, CytO17 should not be included in the linear models that follow.

Analysis of data set:

Question 2: What cytokines in the pre-transplant blood sample specifically correlate with:

  1. Immune cell infiltration after transplant? (Inflammation score)

Since we are looking at cytokines in pre-transplant blood samples, we first filter for pre-operation data.

# filter for preoperation data
dataPO <- dataclean %>% 
  filter(visittype == "PO")

# crosstabs on preoperation data
xtabs(~ Inflammation, data = dataPO)
## Inflammation
##  0  1  2  3 
##  9 59 26  6
xtabs(~ Necrosis, data = dataPO)
## Necrosis
##  0  1  2  3  4 
## 25 37 24 12  2

We reassign the numerical inflammation scores to reflect negative and positive scores.

dataPO1 <- dataPO

# reassigning inflammation scores 0=negative and 2=positive
dataPO1$Inflammation[dataPO$Inflammation == 0] <- 0
dataPO1$Inflammation[dataPO$Inflammation == 1] <- 0

dataPO1$Inflammation[dataPO$Inflammation == 2] <- 2
dataPO1$Inflammation[dataPO$Inflammation == 3] <- 2
dataPO1$Inflammation[dataPO$Inflammation == 4] <- 2


# reassigning necrosis scores 0=negative and 2=positive
dataPO1$Necrosis[dataPO$Necrosis == 0] <- 0
dataPO1$Necrosis[dataPO$Necrosis == 1] <- 0

dataPO1$Necrosis[dataPO$Necrosis == 2] <- 2
dataPO1$Necrosis[dataPO$Necrosis == 3] <- 2
dataPO1$Necrosis[dataPO$Necrosis == 4] <- 2

# crosstabs to check that reassigned values match original data
xtabs(~ Inflammation, data = dataPO1)
## Inflammation
##  0  2 
## 68 32
xtabs(~ Necrosis, data = dataPO1)
## Necrosis
##  0  2 
## 62 38

We attempt an ordinal logistic regression model with the numeric inflammation scores. The confidence interval for cytokine 23 does not include 0, which suggests that it may be statistically significant.

library(MASS)
# releveled inflammation as factor
dataPO$Inflammation <- factor(dataPO$Inflammation, levels = c("0", "1", "2", "3", "4"))
dataPO$Inflammation
##   [1] 1 1 1 1 2 2 2 1 3 1 1 1 3 0 2 3 2 2 1 1 2 1 2 1 2 1 1 1 2 1 2 2 1 1 1 1 1
##  [38] 1 1 1 3 1 2 1 1 1 2 1 3 1 1 2 2 1 1 1 3 1 1 1 2 1 2 1 2 2 1 2 1 1 0 1 1 1
##  [75] 1 1 1 2 1 1 0 1 2 1 1 1 0 2 1 2 1 1 0 2 0 0 0 1 0 1
## Levels: 0 1 2 3 4
# ordinal logistic regression model
m1 <- polr(Inflammation ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO, Hess=TRUE)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(m1)
## Call:
## polr(formula = Inflammation ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + 
##     CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + 
##     CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + 
##     CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO, Hess = TRUE)
## 
## Coefficients:
##             Value Std. Error  t value
## CYTO1   0.0990972  0.0607496  1.63124
## CYTO2  -0.0041909  0.0037136 -1.12854
## CYTO3  -0.0015958  0.0012598 -1.26675
## CYTO4  -0.0185309  0.0103424 -1.79174
## CYTO5   0.0010193  0.0018938  0.53824
## CYTO6   0.0011034  0.0011182  0.98678
## CYTO7   0.0012252  0.0012310  0.99532
## CYTO8   0.0007049  0.0140244  0.05026
## CYTO9   0.0081579  0.0136709  0.59674
## CYTO10 -0.0244423  0.0350566 -0.69722
## CYTO11 -0.0026796  0.0015254 -1.75668
## CYTO12 -0.0264420  0.0223970 -1.18060
## CYTO13 -0.0005690  0.0007607 -0.74808
## CYTO14 -0.0046667  0.0035552 -1.31264
## CYTO15  0.0227108  0.0506707  0.44820
## CYTO16  0.2700851  0.1510910  1.78757
## CYTO18 -0.0003892  0.0002172 -1.79196
## CYTO19  0.0405524  0.0356837  1.13644
## CYTO20 -0.0008111  0.0007264 -1.11657
## CYTO21  0.0233837  0.0277171  0.84366
## CYTO22 -0.0004304  0.0024876 -0.17303
## CYTO23  0.0578270  0.0223748  2.58447
## 
## Intercepts:
##     Value   Std. Error t value
## 0|1 -2.0028  0.6858    -2.9204
## 1|2  1.8115  0.6294     2.8781
## 2|3  4.4719  0.8403     5.3217
## 3|4 23.4397  0.8402    27.8988
## 
## Residual Deviance: 170.8675 
## AIC: 222.8675
## store table
(ctable1 <- coef(summary(m1)))
##                Value   Std. Error     t value
## CYTO1   0.0990972489 0.0607495853  1.63124157
## CYTO2  -0.0041909214 0.0037135846 -1.12853801
## CYTO3  -0.0015957931 0.0012597556 -1.26674817
## CYTO4  -0.0185309221 0.0103423908 -1.79174453
## CYTO5   0.0010193281 0.0018938064  0.53824307
## CYTO6   0.0011033743 0.0011181528  0.98678309
## CYTO7   0.0012252488 0.0012310094  0.99532041
## CYTO8   0.0007049127 0.0140243958  0.05026332
## CYTO9   0.0081579103 0.0136709009  0.59673538
## CYTO10 -0.0244423105 0.0350565760 -0.69722469
## CYTO11 -0.0026795576 0.0015253564 -1.75667638
## CYTO12 -0.0264420415 0.0223970444 -1.18060406
## CYTO13 -0.0005690480 0.0007606737 -0.74808429
## CYTO14 -0.0046667285 0.0035552334 -1.31263633
## CYTO15  0.0227107980 0.0506706904  0.44820384
## CYTO16  0.2700850984 0.1510910309  1.78756540
## CYTO18 -0.0003892376 0.0002172136 -1.79195821
## CYTO19  0.0405524257 0.0356837272  1.13644030
## CYTO20 -0.0008110602 0.0007263833 -1.11657342
## CYTO21  0.0233837014 0.0277170687  0.84365709
## CYTO22 -0.0004304278 0.0024876217 -0.17302785
## CYTO23  0.0578269674 0.0223747943  2.58446923
## 0|1    -2.0028022044 0.6857895393 -2.92043271
## 1|2     1.8114728183 0.6293997691  2.87809578
## 2|3     4.4719234340 0.8403158047  5.32171763
## 3|4    23.4396649994 0.8401688913 27.89875374
## calculate and store p values
p1 <- pnorm(abs(ctable1[, "t value"]), lower.tail = FALSE) * 2

## combined table
(ctable1 <- cbind(ctable1, "p value" = p1))
##                Value   Std. Error     t value       p value
## CYTO1   0.0990972489 0.0607495853  1.63124157  1.028394e-01
## CYTO2  -0.0041909214 0.0037135846 -1.12853801  2.590928e-01
## CYTO3  -0.0015957931 0.0012597556 -1.26674817  2.052453e-01
## CYTO4  -0.0185309221 0.0103423908 -1.79174453  7.317390e-02
## CYTO5   0.0010193281 0.0018938064  0.53824307  5.904093e-01
## CYTO6   0.0011033743 0.0011181528  0.98678309  3.237490e-01
## CYTO7   0.0012252488 0.0012310094  0.99532041  3.195805e-01
## CYTO8   0.0007049127 0.0140243958  0.05026332  9.599126e-01
## CYTO9   0.0081579103 0.0136709009  0.59673538  5.506841e-01
## CYTO10 -0.0244423105 0.0350565760 -0.69722469  4.856622e-01
## CYTO11 -0.0026795576 0.0015253564 -1.75667638  7.897299e-02
## CYTO12 -0.0264420415 0.0223970444 -1.18060406  2.377601e-01
## CYTO13 -0.0005690480 0.0007606737 -0.74808429  4.544093e-01
## CYTO14 -0.0046667285 0.0035552334 -1.31263633  1.893055e-01
## CYTO15  0.0227107980 0.0506706904  0.44820384  6.540061e-01
## CYTO16  0.2700850984 0.1510910309  1.78756540  7.384615e-02
## CYTO18 -0.0003892376 0.0002172136 -1.79195821  7.313966e-02
## CYTO19  0.0405524257 0.0356837272  1.13644030  2.557723e-01
## CYTO20 -0.0008110602 0.0007263833 -1.11657342  2.641768e-01
## CYTO21  0.0233837014 0.0277170687  0.84365709  3.988611e-01
## CYTO22 -0.0004304278 0.0024876217 -0.17302785  8.626295e-01
## CYTO23  0.0578269674 0.0223747943  2.58446923  9.752901e-03
## 0|1    -2.0028022044 0.6857895393 -2.92043271  3.495457e-03
## 1|2     1.8114728183 0.6293997691  2.87809578  4.000836e-03
## 2|3     4.4719234340 0.8403158047  5.32171763  1.027920e-07
## 3|4    23.4396649994 0.8401688913 27.89875374 2.762585e-171
# CIs assuming normality
confint.default(m1)
##                2.5 %       97.5 %
## CYTO1  -0.0199697503 2.181642e-01
## CYTO2  -0.0114694136 3.087571e-03
## CYTO3  -0.0040648686 8.732825e-04
## CYTO4  -0.0388016355 1.739791e-03
## CYTO5  -0.0026924641 4.731120e-03
## CYTO6  -0.0010881650 3.294914e-03
## CYTO7  -0.0011874853 3.637983e-03
## CYTO8  -0.0267823980 2.819222e-02
## CYTO9  -0.0186365631 3.495238e-02
## CYTO10 -0.0931519370 4.426732e-02
## CYTO11 -0.0056692013 3.100861e-04
## CYTO12 -0.0703394419 1.745536e-02
## CYTO13 -0.0020599411 9.218450e-04
## CYTO14 -0.0116348580 2.301401e-03
## CYTO15 -0.0766019302 1.220235e-01
## CYTO16 -0.0260478806 5.662181e-01
## CYTO18 -0.0008149684 3.649313e-05
## CYTO19 -0.0293863945 1.104912e-01
## CYTO20 -0.0022347453 6.126248e-04
## CYTO21 -0.0309407550 7.770816e-02
## CYTO22 -0.0053060769 4.445221e-03
## CYTO23  0.0139731765 1.016808e-01
#Using step function 
step_m1 <- step(m1, trace = F)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

We also attempt a binary response model with the positive/negative inflammation scores, which seems to be a better fitted model. The confidence intervals for cytokine 4, 9, 10, and 23 do not include 0, which suggests that it may be statistically significant.

xtabs(~ Inflammation, data = dataclean)
## Inflammation
##   0   1   2   3 
##  27 177  78  18
xtabs(~ Necrosis, data = dataclean)
## Necrosis
##   0   1   2   3   4 
##  75 111  72  36   6
# split outcomes (Inflammation or Necrosis Score) with 0-1 vs 2+ as -/+
# try two types of models: ordered factor outcome and binary outcome (try both)
# releveled inflammation as factor
dataPO1$Inflammation <- factor(dataPO1$Inflammation)

# binary response model
m3 <- glm(Inflammation ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO1, family = "binomial")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(m3)
## 
## Call:
## glm(formula = Inflammation ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + 
##     CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + 
##     CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + 
##     CYTO20 + CYTO21 + CYTO22 + CYTO23, family = "binomial", data = dataPO1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3855  -0.7277  -0.4042   0.5838   2.5194  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -1.9594398  0.7981657  -2.455  0.01409 * 
## CYTO1        0.1654207  0.0899441   1.839  0.06589 . 
## CYTO2       -0.0012708  0.0051856  -0.245  0.80641   
## CYTO3       -0.0015651  0.0014822  -1.056  0.29101   
## CYTO4       -0.0330060  0.0158206  -2.086  0.03695 * 
## CYTO5        0.0069897  0.0120757   0.579  0.56271   
## CYTO6        0.0010706  0.0013156   0.814  0.41575   
## CYTO7        0.0010141  0.0012318   0.823  0.41033   
## CYTO8        0.0313095  0.0292942   1.069  0.28516   
## CYTO9        0.0409934  0.0207831   1.972  0.04856 * 
## CYTO10      -0.1728190  0.0837291  -2.064  0.03902 * 
## CYTO11      -0.0043394  0.0026628  -1.630  0.10317   
## CYTO12      -0.0474168  0.0549345  -0.863  0.38805   
## CYTO13      -0.0008085  0.0007021  -1.152  0.24951   
## CYTO14      -0.0067992  0.0054303  -1.252  0.21054   
## CYTO15      -0.1511398  0.1841940  -0.821  0.41190   
## CYTO16       0.5637371  0.3519016   1.602  0.10916   
## CYTO18      -0.0012869  0.0010949  -1.175  0.23985   
## CYTO19       0.0608659  0.0488801   1.245  0.21306   
## CYTO20      -0.0020995  0.0025164  -0.834  0.40410   
## CYTO21       0.0050385  0.0412417   0.122  0.90276   
## CYTO22       0.0039077  0.0033793   1.156  0.24754   
## CYTO23       0.0879739  0.0311802   2.821  0.00478 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 125.374  on 99  degrees of freedom
## Residual deviance:  87.069  on 77  degrees of freedom
## AIC: 133.07
## 
## Number of Fisher Scoring iterations: 8
# CIs assuming normality
confint.default(m3)
##                     2.5 %        97.5 %
## (Intercept) -3.5238158763 -0.3950636782
## CYTO1       -0.0108665844  0.3417079653
## CYTO2       -0.0114343560  0.0088927965
## CYTO3       -0.0044702025  0.0013400241
## CYTO4       -0.0640138551 -0.0019981919
## CYTO5       -0.0166781949  0.0306575063
## CYTO6       -0.0015078163  0.0036490653
## CYTO7       -0.0014000859  0.0034283666
## CYTO8       -0.0261060666  0.0887251165
## CYTO9        0.0002592515  0.0817276223
## CYTO10      -0.3369249713 -0.0087129802
## CYTO11      -0.0095583507  0.0008795283
## CYTO12      -0.1550864816  0.0602529488
## CYTO13      -0.0021845399  0.0005675757
## CYTO14      -0.0174424859  0.0038440589
## CYTO15      -0.5121534497  0.2098738316
## CYTO16      -0.1259773237  1.2534514546
## CYTO18      -0.0034329049  0.0008590563
## CYTO19      -0.0349374532  0.1566691788
## CYTO20      -0.0070314967  0.0028325405
## CYTO21      -0.0757937711  0.0858707302
## CYTO22      -0.0027156585  0.0105310998
## CYTO23       0.0268619062  0.1490859247

We also try a stepwise regression model, which gives the best fittedd model. Cytokines 9, 16, 18, and 23 have p-values less than 0.05, which suggests that they may specifically correlate with immune cell infiltration after transplant.

library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
library(leaps)
library(MASS)

# Stepwise regression model
step.model <- stepAIC(m3, direction = "both", trace = FALSE)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
# step.model <- full.model %>% stepAIC(trace = FALSE)

summary(step.model)
## 
## Call:
## glm(formula = Inflammation ~ CYTO4 + CYTO6 + CYTO9 + CYTO10 + 
##     CYTO11 + CYTO16 + CYTO18 + CYTO23, family = "binomial", data = dataPO1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5199  -0.8133  -0.5408   0.8247   2.1625  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.5454002  0.6126235  -2.523   0.0116 *
## CYTO4       -0.0167166  0.0113717  -1.470   0.1416  
## CYTO6        0.0014866  0.0008917   1.667   0.0955 .
## CYTO9        0.0406465  0.0157946   2.573   0.0101 *
## CYTO10      -0.0771683  0.0407839  -1.892   0.0585 .
## CYTO11      -0.0033005  0.0017340  -1.903   0.0570 .
## CYTO16       0.3184029  0.1309310   2.432   0.0150 *
## CYTO18      -0.0013273  0.0005463  -2.429   0.0151 *
## CYTO23       0.0551877  0.0226520   2.436   0.0148 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 125.37  on 99  degrees of freedom
## Residual deviance: 101.14  on 91  degrees of freedom
## AIC: 119.14
## 
## Number of Fisher Scoring iterations: 7
  1. Dead cells in the liver after transplant? (Necrosis score)

We attempt an ordinal logistic regression model on the numerical necrosis scores. All of the confidence intervals include 0, so none of the cytokines appear to be statistically significant.

library(MASS)
# releveled necrosis as factor
dataPO$Necrosis <- factor(dataPO$Necrosis, levels = c("0", "1", "2", "3", "4"))
dataPO$Necrosis
##   [1] 0 0 0 0 1 4 2 1 3 1 2 1 1 0 1 4 3 2 3 1 3 2 2 0 1 0 1 2 2 2 2 2 1 3 2 1 1
##  [38] 1 1 1 2 1 2 1 1 1 2 3 2 1 0 3 2 1 2 1 3 2 2 0 0 1 1 3 3 2 0 1 1 1 0 0 3 0
##  [75] 3 0 0 2 2 0 0 1 1 2 0 0 0 1 1 1 1 1 0 1 0 1 0 0 2 1
## Levels: 0 1 2 3 4
# ordinal logistic regression model
m2 <- polr(Necrosis ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5  + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO, Hess=TRUE)
summary(m2)
## Warning in sqrt(diag(vc)): NaNs produced
## Call:
## polr(formula = Necrosis ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + 
##     CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + 
##     CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + CYTO20 + 
##     CYTO21 + CYTO22 + CYTO23, data = dataPO, Hess = TRUE)
## 
## Coefficients:
##             Value Std. Error  t value
## CYTO1   3.192e-02  0.0565505  0.56446
## CYTO2  -2.545e-03        NaN      NaN
## CYTO3   1.250e-03  0.0007780  1.60670
## CYTO4  -1.025e-02  0.0098028 -1.04552
## CYTO5  -1.687e-03  0.0016441 -1.02605
## CYTO6  -1.043e-03  0.0009796 -1.06516
## CYTO7   1.455e-03  0.0012184  1.19419
## CYTO8   4.492e-03        NaN      NaN
## CYTO9  -2.141e-03        NaN      NaN
## CYTO10 -1.346e-02        NaN      NaN
## CYTO11  1.323e-03  0.0014214  0.93054
## CYTO12  6.778e-03  0.0146035  0.46417
## CYTO13 -6.015e-04  0.0011060 -0.54383
## CYTO14 -2.169e-03  0.0035172 -0.61663
## CYTO15 -3.605e-02        NaN      NaN
## CYTO16  7.691e-02        NaN      NaN
## CYTO18 -4.726e-04        NaN      NaN
## CYTO19  2.952e-02        NaN      NaN
## CYTO20 -1.874e-04        NaN      NaN
## CYTO21  3.212e-02        NaN      NaN
## CYTO22  4.788e-05  0.0017506  0.02735
## CYTO23  2.922e-02  0.0192031  1.52143
## 
## Intercepts:
##     Value   Std. Error t value
## 0|1 -0.5193  0.5504    -0.9434
## 1|2  1.2675  0.5629     2.2519
## 2|3  2.7644  0.6426     4.3018
## 3|4  4.9783  0.9831     5.0637
## 
## Residual Deviance: 259.1091 
## AIC: 311.1091
## store table
(ctable2 <- coef(summary(m2)))
## Warning in sqrt(diag(vc)): NaNs produced
##                Value   Std. Error     t value
## CYTO1   3.192063e-02 0.0565505255  0.56446210
## CYTO2  -2.544541e-03          NaN         NaN
## CYTO3   1.249995e-03 0.0007779874  1.60670338
## CYTO4  -1.024897e-02 0.0098027847 -1.04551616
## CYTO5  -1.686893e-03 0.0016440617 -1.02605191
## CYTO6  -1.043408e-03 0.0009795748 -1.06516424
## CYTO7   1.454960e-03 0.0012183636  1.19419169
## CYTO8   4.492127e-03          NaN         NaN
## CYTO9  -2.140751e-03          NaN         NaN
## CYTO10 -1.345868e-02          NaN         NaN
## CYTO11  1.322672e-03 0.0014214048  0.93053825
## CYTO12  6.778477e-03 0.0146034506  0.46416955
## CYTO13 -6.014651e-04 0.0011059715 -0.54383420
## CYTO14 -2.168817e-03 0.0035172148 -0.61662903
## CYTO15 -3.604522e-02          NaN         NaN
## CYTO16  7.691065e-02          NaN         NaN
## CYTO18 -4.725605e-04          NaN         NaN
## CYTO19  2.952395e-02          NaN         NaN
## CYTO20 -1.874131e-04          NaN         NaN
## CYTO21  3.211881e-02          NaN         NaN
## CYTO22  4.788159e-05 0.0017505990  0.02735155
## CYTO23  2.921612e-02 0.0192031158  1.52142586
## 0|1    -5.192508e-01 0.5504176604 -0.94337602
## 1|2     1.267495e+00 0.5628545961  2.25190446
## 2|3     2.764412e+00 0.6426239226  4.30175733
## 3|4     4.978260e+00 0.9831229354  5.06372054
## calculate and store p values
p2 <- pnorm(abs(ctable2[, "t value"]), lower.tail = FALSE) * 2

## combined table
(ctable2 <- cbind(ctable2, "p value" = p2))
##                Value   Std. Error     t value      p value
## CYTO1   3.192063e-02 0.0565505255  0.56446210 5.724397e-01
## CYTO2  -2.544541e-03          NaN         NaN          NaN
## CYTO3   1.249995e-03 0.0007779874  1.60670338 1.081195e-01
## CYTO4  -1.024897e-02 0.0098027847 -1.04551616 2.957845e-01
## CYTO5  -1.686893e-03 0.0016440617 -1.02605191 3.048671e-01
## CYTO6  -1.043408e-03 0.0009795748 -1.06516424 2.868016e-01
## CYTO7   1.454960e-03 0.0012183636  1.19419169 2.324030e-01
## CYTO8   4.492127e-03          NaN         NaN          NaN
## CYTO9  -2.140751e-03          NaN         NaN          NaN
## CYTO10 -1.345868e-02          NaN         NaN          NaN
## CYTO11  1.322672e-03 0.0014214048  0.93053825 3.520925e-01
## CYTO12  6.778477e-03 0.0146034506  0.46416955 6.425263e-01
## CYTO13 -6.014651e-04 0.0011059715 -0.54383420 5.865556e-01
## CYTO14 -2.168817e-03 0.0035172148 -0.61662903 5.374794e-01
## CYTO15 -3.604522e-02          NaN         NaN          NaN
## CYTO16  7.691065e-02          NaN         NaN          NaN
## CYTO18 -4.725605e-04          NaN         NaN          NaN
## CYTO19  2.952395e-02          NaN         NaN          NaN
## CYTO20 -1.874131e-04          NaN         NaN          NaN
## CYTO21  3.211881e-02          NaN         NaN          NaN
## CYTO22  4.788159e-05 0.0017505990  0.02735155 9.781793e-01
## CYTO23  2.921612e-02 0.0192031158  1.52142586 1.281530e-01
## 0|1    -5.192508e-01 0.5504176604 -0.94337602 3.454886e-01
## 1|2     1.267495e+00 0.5628545961  2.25190446 2.432831e-02
## 2|3     2.764412e+00 0.6426239226  4.30175733 1.694488e-05
## 3|4     4.978260e+00 0.9831229354  5.06372054 4.111522e-07
# CIs assuming normality
confint.default(m2)
## Warning in sqrt(diag(vcov(object))): NaNs produced
##                2.5 %       97.5 %
## CYTO1  -0.0789163649 0.1427576217
## CYTO2            NaN          NaN
## CYTO3  -0.0002748323 0.0027748222
## CYTO4  -0.0294620747 0.0089641351
## CYTO5  -0.0049091943 0.0015354091
## CYTO6  -0.0029633395 0.0008765233
## CYTO7  -0.0009329891 0.0038429084
## CYTO8            NaN          NaN
## CYTO9            NaN          NaN
## CYTO10           NaN          NaN
## CYTO11 -0.0014632307 0.0041085737
## CYTO12 -0.0218437601 0.0354007144
## CYTO13 -0.0027691295 0.0015661992
## CYTO14 -0.0090624312 0.0047247976
## CYTO15           NaN          NaN
## CYTO16           NaN          NaN
## CYTO18           NaN          NaN
## CYTO19           NaN          NaN
## CYTO20           NaN          NaN
## CYTO21           NaN          NaN
## CYTO22 -0.0033832294 0.0034789926
## CYTO23 -0.0084212983 0.0668535326

We also attempt a binary response model on the positive/negative necrosis scores, which gives a better fitted model. Again, all of the confidence intervals include 0, so none of the cytokines appear to be statistically significant.

# releveled necrosis as factor
dataPO1$Necrosis <- factor(dataPO1$Necrosis)

#binary response model
m4 <- glm(Necrosis ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO17 + CYTO18 + CYTO19 + CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO1, family = "binomial")
summary(m4)
## 
## Call:
## glm(formula = Necrosis ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + 
##     CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + 
##     CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO17 + CYTO18 + CYTO19 + 
##     CYTO20 + CYTO21 + CYTO22 + CYTO23, family = "binomial", data = dataPO1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4905  -0.9435  -0.5007   1.0971   2.1463  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -0.0304548  1.2074987  -0.025   0.9799  
## CYTO1        0.1051974  0.0757430   1.389   0.1649  
## CYTO2       -0.0072987  0.0041645  -1.753   0.0797 .
## CYTO3        0.0009495  0.0013069   0.727   0.4675  
## CYTO4        0.0011742  0.0107202   0.110   0.9128  
## CYTO5       -0.0124740  0.0208794  -0.597   0.5502  
## CYTO6       -0.0014713  0.0013816  -1.065   0.2869  
## CYTO7        0.0017876  0.0011727   1.524   0.1274  
## CYTO8       -0.0144491  0.0245025  -0.590   0.5554  
## CYTO9        0.0097545  0.0163221   0.598   0.5501  
## CYTO10       0.0150727  0.0483409   0.312   0.7552  
## CYTO11      -0.0001271  0.0018280  -0.070   0.9446  
## CYTO12      -0.0138525  0.0240642  -0.576   0.5649  
## CYTO13      -0.0008963  0.0006146  -1.458   0.1448  
## CYTO14      -0.0085112  0.0047045  -1.809   0.0704 .
## CYTO15       0.0320340  0.0839708   0.381   0.7028  
## CYTO16       0.1078550  0.1292101   0.835   0.4039  
## CYTO17      -1.0104245  0.7120463  -1.419   0.1559  
## CYTO18      -0.0001255  0.0002660  -0.472   0.6372  
## CYTO19       0.0400669  0.0378243   1.059   0.2895  
## CYTO20      -0.0011347  0.0009620  -1.180   0.2382  
## CYTO21       0.0379885  0.0344121   1.104   0.2696  
## CYTO22      -0.0007377  0.0027145  -0.272   0.7858  
## CYTO23       0.0322030  0.0252215   1.277   0.2017  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 132.81  on 99  degrees of freedom
## Residual deviance: 110.52  on 76  degrees of freedom
## AIC: 158.52
## 
## Number of Fisher Scoring iterations: 7
# CIs assuming normality
confint.default(m4)
##                     2.5 %       97.5 %
## (Intercept) -2.3971086348 2.3361991026
## CYTO1       -0.0432560885 0.2536508289
## CYTO2       -0.0154610355 0.0008635451
## CYTO3       -0.0016118691 0.0035109007
## CYTO4       -0.0198369429 0.0221854280
## CYTO5       -0.0533967639 0.0284488620
## CYTO6       -0.0041791609 0.0012364694
## CYTO7       -0.0005108973 0.0040861837
## CYTO8       -0.0624731722 0.0335750531
## CYTO9       -0.0222361432 0.0417451927
## CYTO10      -0.0796737313 0.1098191153
## CYTO11      -0.0037099215 0.0034557073
## CYTO12      -0.0610173986 0.0333124625
## CYTO13      -0.0021010175 0.0003083590
## CYTO14      -0.0177319133 0.0007095249
## CYTO15      -0.1325457134 0.1966137291
## CYTO16      -0.1453920737 0.3611021632
## CYTO17      -2.4060095573 0.3851604671
## CYTO18      -0.0006468930 0.0003959766
## CYTO19      -0.0340673446 0.1142011674
## CYTO20      -0.0030202375 0.0007507667
## CYTO21      -0.0294579473 0.1054350298
## CYTO22      -0.0060580481 0.0045825583
## CYTO23      -0.0172302817 0.0816363534

We also try a stepwise regression model, which gives the best fitted model. Cytokine 1 has a p-value less than 0.05, which suggests that this cytokine specifically correlates with the dead cells in the liver after transplant.

library(caret)
library(leaps)
library(MASS)

# Stepwise regression model
step.model <- stepAIC(m4, direction = "both", 
                      trace = FALSE)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(step.model)
## 
## Call:
## glm(formula = Necrosis ~ CYTO1 + CYTO2 + CYTO13 + CYTO17 + CYTO18 + 
##     CYTO23, family = "binomial", data = dataPO1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5453  -0.9695  -0.6564   1.2369   1.8411  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -0.1873922  0.9237729  -0.203   0.8392  
## CYTO1        0.0902847  0.0455846   1.981   0.0476 *
## CYTO2       -0.0041882  0.0024074  -1.740   0.0819 .
## CYTO13      -0.0005261  0.0004327  -1.216   0.2240  
## CYTO17      -0.8517493  0.5849126  -1.456   0.1453  
## CYTO18      -0.0003877  0.0003022  -1.283   0.1996  
## CYTO23       0.0350102  0.0193580   1.809   0.0705 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 132.81  on 99  degrees of freedom
## Residual deviance: 119.89  on 93  degrees of freedom
## AIC: 133.89
## 
## Number of Fisher Scoring iterations: 6

Interpretation of results in terms of the original research hypotheses: